Wstęp

Celem poniższego projektu jest analiza zbioru danych dotyczących pomiarów śledzi atlantyckich i warunków w jakich żyją, zbieranych w ciągu ostatnich 60 lat, i próba znalezienia przyczyn stopniowego spadku ich długości w ostatnich latach. W procesie tej analizy zbadano ewentualne zależności pomiędzy atrybutami zbioru i utworzono klasyfikator dokonujący predykcji rozmiaru śledzi na podstawie pozostałych pomiarów. W wyniku poniższych kroków sformułowano trzy główne czynniki wpływające na spadek długości śledzi w ostatnich latach.

Przygotowanie środowiska

Poniższy kod ładuje wymagane biblioteki i zapewnia powtarzalność wyników:

library(ggplot2)
library(plotly)
library(VIM)
library(knitr)
library(ggcorrplot)
library(caret)       
library(randomForest)
set.seed(42)

Zbiór danych

Zbiór danych składa się z pomiarów śledzi i warunków w jakich żyją, zebranych przez ostatnich 60 lat. Dane były pobierane z połowów komercyjnych jednostek. W ramach połowu jednej jednostki losowo wybierano od 50 do 100 sztuk trzyletnich śledzi.

Kolejne kolumny w zbiorze danych to:

Oryginalny zbiór

Wczytywanie oryginalnego zbioru danych z wyspecyfikowaniem typu danych w kolumnach:

sample <- read.csv("sledzie.csv", nrows = 100)
df_classes <- sapply(sample, class) 
df_classes <- replace(df_classes, df_classes=="factor", "numeric")
all_df <- read.csv("sledzie.csv", colClasses = df_classes, fill=TRUE, na.string=c("NA", "?"))

Próbka danych ze zbioru:

X length cfin1 cfin2 chel1 chel2 lcop1 lcop2 fbar recr cumf totaln sst sal xmonth nao
0 23.0 0.02778 0.27785 2.46875 NA 2.54787 26.35881 0.356 482831 0.3059879 267380.8 14.30693 35.51234 7 2.8
1 22.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831 0.3059879 267380.8 14.30693 35.51234 7 2.8
2 25.0 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831 0.3059879 267380.8 14.30693 35.51234 7 2.8
3 25.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831 0.3059879 267380.8 14.30693 35.51234 7 2.8
4 24.0 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831 0.3059879 267380.8 14.30693 35.51234 7 2.8
5 22.0 0.02778 0.27785 2.46875 21.43548 2.54787 NA 0.356 482831 0.3059879 267380.8 14.30693 35.51234 7 2.8
6 24.0 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831 0.3059879 267380.8 14.30693 35.51234 7 2.8
7 23.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831 0.3059879 267380.8 14.30693 35.51234 7 2.8
8 22.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831 0.3059879 267380.8 14.30693 35.51234 7 2.8
9 22.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831 0.3059879 267380.8 14.30693 35.51234 7 2.8

Uzupełniony zbiór danych

Oryginalny zbiór składa się z 52582 pomiarów, z czego 10094 (około 19%) zawiera brakujące wartości:

##   n_samples n_missing_values percent_missing_values
## 1     52582            10094              0.1919668

Poniżej zostały zaprezentowane szczegółowe informacje na temat brakujących wartości pomiarów:

## 
##  Variables sorted by number of missings: 
##  Variable      Count
##     lcop1 0.03143661
##     lcop2 0.03025750
##       sst 0.03012438
##     cfin1 0.03006732
##     chel2 0.02959188
##     chel1 0.02957286
##     cfin2 0.02921152
##         X 0.00000000
##    length 0.00000000
##      fbar 0.00000000
##      recr 0.00000000
##      cumf 0.00000000
##    totaln 0.00000000
##       sal 0.00000000
##    xmonth 0.00000000
##       nao 0.00000000

Można zauważyć, że braki w pomiarach dotyczą 7 atrybutów - dostępności poszczególnych gatunków planktonu oraz temperatury przy powierzchni wody.

Możliwym rozwiązaniem tego problemu byłoby użycie metody k-nn, aby uzupełniać brakujące dane w pomiarze na podstawie tych znajdujących się w najpodobniejszych do niego pomiarach.

Po dokładniejszym przyjrzeniu się danym, okazuje się, że w przypadku większości takich wierszy można łatwo uzupełnić te braki, ponieważ otaczające je pomiary, pochodzące z tego samego połowu, są bardzo podobne lub identyczne, mogą więc być wykorzystane do wypełnienia braków, bez uszczerbku na ogólnej jakości zbioru.

Ponieważ dane są uporządkowane chronologicznie można założyć, że jeśli dwa sąsiednie pomiary mają taką samą wartość atrybutu xmonth to mogą posłużyć do uzupełnienia braków sąsiada. Aby uprościć proces, w przypadku, jeśli sąsiednie pomiary pochodzą już z innego miesiąca lub mają takie same braki jak pomiar aktualnie uzupełniany to wartość brakująca jest nadpisywana przez średnią dla danego atrybutu w całym zbiorze.

Brakujące dane uzupełnia następujący kod:

missing_vals <- apply(all_df, 1, anyNA)

if(missing_vals[1]){
  for (x in 1:ncol(all_df)) {
    if (is.na(all_df[1, x])) {
      all_df[1, x] <- all_df[2, x]
    }
  }
}

for(i in 2:nrow(all_df)-1){
  if(missing_vals[i]){
    for (x in 1:ncol(all_df)) {
      if (is.na(all_df[i, x])) {
        if(!is.na(all_df[i - 1, x]) && (all_df$xmonth[i] == all_df$xmonth[i - 1])){
          all_df[i, x] <- all_df[i - 1, x]
        }else if(!is.na(all_df[i + 1, x]) && (all_df$xmonth[i] == all_df$xmonth[i + 1])){
          all_df[i, x] <- all_df[i + 1, x]
        }else{
          all_df[i, x] <- mean(all_df[,x], na.rm = TRUE)
        }
      }
    }
  }
}

if(missing_vals[nrow(all_df)]){
  for (x in 1:ncol(all_df)) {
    if (is.na(all_df[nrow(all_df), x])) {
      all_df[nrow(all_df), x] <- all_df[nrow(all_df)-1, x]
    }
  }
}

Niestety atrybut xmonth nie wystarczy by jednoznacznie oddzielić kolejne lata pomiarów, ponieważ jego wartość zmienia się częściej niż by to wynikało z rzeczywistego okresu zbierania pomiarów:

months <- 0
curr_month <-0
for(i in 1:nrow(all_df)){
  if(curr_month != all_df$xmonth[i]){
    months <- months + 1
    curr_month <- all_df$xmonth[i]
  }
}
months/12
## [1] 317.5833

W związku z tym zostały utworzone sztuczne atrybuty year oraz years_group, dzielące zbiór odpowiednio na 60 i 10 równej wielkości grup. Jest to podział dokonany na potrzeby analizy i wizualizacji danych.

Analiza wartości atrybutów

Długość wyławianych śledzi

Przebieg poniższego wykresu obrazuje stopniowy spadek rozmiaru śledzia z czasem:

Średnia dostępność planktonu

Na wykresach można zaobserwować zmniejszenie się dostępności wszystkich gatunków planktonu w ostatnich latach, co pokrywa się z mniejszymi rozmiarami osiąganymi przez śledzie. Zależność wydaje się być szczególnie widoczna w przypadku chel1 i lcop1. Potwierdza to intuicyjny wniosek, że dostępność pokarmu może mieć znaczący wpływ na rozmiary ryb.

Dane dotyczące połowu śledzi

Powyższe wykresy nie obrazują wyraźnych zależności. Szczególnie mała instensywność połowów przypada zarówno na okres osiągania przez śledzie największych jak i najmniejszych odnotowanych rozmiarów. Dodatkowo ze względu na problem z wyraźnym podziałem pomiarów na lata w zbiorze danych atrybuty “roczne” stają się mniej wiarygodne.

Dane dotyczące oceanu

Chociaż poziom zasolenia nie wydaje się mieć wielkiego wpływu na rozmiary śledzi, to wykresy sugerują znaczący wpływ temperatury wody i oscylacji północnoatlantyckiej w tej kwestii.

Badanie korelacji

Badanie wykazało istnienie korelacji między długością śledzi a dostępnością niektórych gatunków planktonu (chel1 i lcop1) oraz ułamkiem pozostawionego narybku (fbar), a także odwrotnej korelacji z temperaturą przy powierzchni wody (sst) oraz oscylacją północnoatlantycką (nao). Można więc wnioskować, że na zmniejszenie się rozmiarów osiąganych przez śledzie wpłynęły zmniejszona dostępność pokarmu, niekorzystne warunki środowiska oraz zwiększone natężenie połowów. Dodatkowo można także zauważyć zależności między oscylacją a wspomnianymi miarami dostępności planktonu oraz temperaturą wody.

Klasyfikator

W celu predykcji długości wyłowionej ryby została wykorzystana metoda RandomForest, a przyjętymi miarami oceny jakości klasyfikacji \(R^2\) oraz \(RMSE\). Poniżej zaprezentowany został kod tworzący i uczący model oraz wyniki jego działania.

rf_df <- all_df[,c(2:14,16)]

inTraining <- createDataPartition(y = rf_df$length, p = 0.8, list = FALSE)

training <- rf_df[inTraining, ]
testing <- rf_df[-inTraining, ]

Mtries <- tuneRF(training, training$length, ntree = 30, plot = FALSE)
## mtry = 4  OOB error = 0.07235467 
## Searching left ...
## mtry = 2     OOB error = 0.2741899 
## -2.789527 0.05 
## Searching right ...
## mtry = 8     OOB error = 0.003078526 
## 0.9574523 0.05 
## mtry = 14    OOB error = 0.0002513633 
## 0.9183495 0.05
rfGrid <- expand.grid(mtry = Mtries)
gridCtrl <- trainControl(method = "repeatedcv", number = 2, repeats = 5)

fit <- train(length ~ ., data = training, method = "rf", metric = "RMSE", preProc = c("center", "scale"), trControl = gridCtrl, tuneGrid = rfGrid, ntree = 30, importance=TRUE)

rfClasses <- predict(fit, newdata = testing)
postResample(pred = rfClasses, obs = testing$length)
##      RMSE  Rsquared       MAE 
## 1.1937860 0.4712525 0.9453771

Wykres poniżej wskazuje jako najważniejsze atrybuty związane z dostępnością planktonu, temperaturę blisko powierzchni wody oraz łączną liczbę ryb wyłowionych w ramach połowu:

Wnioski

Na podstawie przeprowadzonej analizy można wyróżnić trzy czynniki, mające znaczny wpływ na spadającą długość śledzi w ostatnich latach:

Warto jednak podkreślić, że śledź atlantycki jest częścią większego ekosystemu i na jego rozwój wpływać będą liczne czynniki, które nie zostały uwzględnione w zbiorze danych. Być może na podstawie szerszego zbioru atrybutów udałoby się znaleźć więcej takich zależności.